#####
# Replication for: "Can Political Speech Foster Tolerance of Immigrants?" by Schleiter, Tavits, and Ward.
# Figure S.7
#####

library(here)
library(data.table)
library(ggplot2)

# source the custom functions
source("functions.R")

# load the pooled data if not already in workspace
if(!exists("pooled")){
  pooled <- fread(here("data", "pooled.csv"))
}

# Create the trichotomous news indicator
pooled[ , news3 := factor(
  fcase(
    news2 < 6, "Low",
    news2 == 6, "Medium",
    news2 == 7, "High"
  ), levels = c("Low", "Medium", "High")
)]


# Helper functions --------------------------------------------------------

# These three functions do the starch of the analysis. The first fits all 9 models, the second extracts the 9 marginal effects per model. The third does the battery of p-values
allCovNewsfits <- 
  function(y, samp, data = pooled, chatty = T){
    if(samp != "Replication"){
      form <- as.formula(paste0(y, "~ CH + NM + CS + CH*news3 + NM*news3 + CS*news3 + female + age + I(age^2) + race2 + pid2 + edu2 + region2 + factor(wave)"))
    } else {
      form <- as.formula(paste0(y, "~ CH + NM + CS + CH*news3 + NM*news3 + CS*news3 + female + age + I(age^2) + race2 + pid2 + edu2 + region2"))
    }
    
    sample_use <- c("Main" = "replication == 0", "Replication" = "replication == 1", "Pooled" = "")[samp]
    
    if(samp != "Pooled"){
      mod <- lm(formula = form, data = data[eval(parse(text = sample_use))])
    } else {
      mod <- lm(formula = form, data = data)
    }
    
    if(chatty) print(summary(mod))
    return(mod)
  }

news_ME_extract <- 
  function(mod, level = .95){
    out <- expand.grid("treatment" = c("CH", "NM", "CS"), "news" = c("L", "M", "H"), stringsAsFactors = F); setDT(out)
    
    out[news == "L" , c("est", "lower", "upper") := .(
      coef(mod)[c("CH", "NM", "CS")],
      confint(mod, level = level)[c("CH", "NM", "CS"), 1],
      confint(mod, level = level)[c("CH", "NM", "CS"), 2]
    )]
    
    for(n in 4:9){
      vars <- c(out$treatment[n], paste0(out$treatment[n], ":news3", ifelse(out$news[n] == "M", "Medium", "High")))
      est <- sum(coef(mod)[vars])
      se <- sqrt(vcov(mod)[vars[1], vars[1]] + vcov(mod)[vars[2], vars[2]] + 2*vcov(mod)[vars[1], vars[2]])
      set(out, i = n, j = "est", value = est)
      set(out, i = n, j = "lower", value = est - 1.96*se)
      set(out, i = n, j = "upper", value = est + 1.96*se)
    }
    
    return(out)
  }

# Fitting -----------------------------------------------------------------


news_fits <- expand.grid("y" = c("immPCA", "imm_neighbors2", "imm_increase2"), "samp" = c("Main", "Replication", "Pooled")); setDT(news_fits)

news_fits[ , fit := Map(allCovNewsfits, y, samp, chatty = F)]
news_me<- news_fits[ , news_ME_extract(mod = fit[[1]]), by = .(y, samp)]


# Plot Prep ---------------------------------------------------------------



news_me[ , yval := c("CH" = 1L, "CS" = 2L, "NM" = 3L)[treatment]]
news_me[ , news2 := factor(
  fcase(
    news == "L", "Low",
    news == "H", "High",
    news == "M", "Medium"
  ),
  levels = c("Low", "Medium", "High")
)]
news_me[ , y2 := factor(y, labels = c("Imm. Index (SD = 1)", "Imm. Neighbors (1-7)", "Increase Imm. (1-7)"))]
news_me[ , samp2 := factor(samp, labels = c("Main Sample", "Replication Sample", "Pooled Sample"))]


# Plotting ----------------------------------------------------------------


(news_plot <- ggplot(news_me, aes(x = est, xmax = upper, xmin = lower, y = yval, color = news2, shape = news2)) +
    geom_vline(xintercept = 0, linetype = "dashed") +
    ggstance::geom_pointrangeh(position = ggstance::position_dodge2v(.7, reverse = T)) +
    facet_grid(y2 ~ samp2) +
    scale_y_continuous(breaks = c(1,2,3), labels = c("Common\nHumanity", "Countering\nStereotypes", "Norms"), limits = c(.5, 3.5), expand = expansion(0,0)) +
    ylab(NULL) +
    xlab("Estimated Treatment Effect") +
    scale_color_grey(name = "News Attention", end = 0.7) +
    scale_shape(name = "News Attention") +
    theme_bw() +
    theme(legend.position = "bottom",
          axis.title.x = element_text(size = 10),
          axis.text = element_text(size = 10),
          strip.text = element_text(size = 10),
          legend.text = element_text(size = 10),
          legend.title = element_text(size = 10))
)
SI_width = 6.50127  

ggsave(here("SI", "figure_s7.pdf"), news_plot, width = SI_width, height = SI_width*1.25)
